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We perform direct numerical simulations of an unstably stratified turbulent channel flow to ad¬ 
dress the effects of buoyancy on the boundary layer dynamics and mean field quantities. We 
systematically span a range of parameters in the space of friction Reynolds number (Rer) and 
Rayleigh number (Ra). Our focus is on deviations from the logarithmic law of the wall due 
to buoyant motion. The effects of convection in the relevant ranges are discussed providing 
measurements of mean profiles of velocity, temperature and Reynolds stresses as well as of the 
friction coefficient. A phenomenological model is proposed and shown to capture the observed 
deviations of the velocity profile in the log-law region from the non-convective case. 


1. Introduction 

Wall bounded turbulence is highly relevant to a variety of engineering systems and naturally 
occurring flows. As such, the importance of the fundamental understanding of the generation of 
turbulent fluctuations by fluid-solid interactions in the boundary layers and their transport and 
interaction with the bulk mean flow ( Pope II200II) is evident. In many geophysical and industrial 
flows, shear flows and wall turbulence are subjected to either stable or unstable thermal stratifica¬ 
tions. Stable stratification is known to inhibit the bursting phenomenon, i.e. depleting the trans¬ 
port of turbulent fluctuations from the wall to the bulk, thus reducing the turbulent drag and in¬ 
creasing the mean flow vel ocity. Many studies have been devoted to this situation when the s trat- 
ifled scalar held i s passive (Ijohansson & Wikstrom 1999^, Pa 2 avassiliou&_Hanra^ 19971) and 


active dArmenio & Sarkar 20021: Gerz. Schumann & Elgobashi 1 1989 : lida. Kasagi & Nagano I 
2002 : GarcIa;Wllalb^^^d^lamo]j20]jJ) and even in presence of non-Oberbeck-Boussinesq 
effects (Izonta. Onorato & Soldati 1 2012h. 

The unstable configuration is also relevant in a variety of instances, such as the physics of the 
atmospheric surface layer superposed to an over-heat ed ground (as, e .g., in summer days, de¬ 
termining the generation of thermo-convective storms ( Bluestein1l2013h ): when buoyancy forces 
are strong enough, thermal convection can occur, altering significantly the unstable boundary 
layer dynamics. Early th eoretical results on turbulent boundary layers under unstable thermal 
stratification date back to ftun^ ( 1932h who pioneered a mixing-length based approach which 
insp ired the later work of Obukhov ( 1946h anticipating the Monin-Obukhov similarity the- 
or v dMonin & Obukhov 1954I). and subsequently revised and compared with experimental data 
by Kader & Yaglomh 199(lh . Lumlev. Zeman & Siess ( 1978h proposed an eddy-damped quasi- 
Gaussian closure able to predict the inversion regions of heat flux profiles in strongly buoyant 
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sheared boundary layers. Numerical results ( lida & K ap^i ||l 997h and expe rimental flume mea¬ 
surements (combined with a spectral equation model) (IKomori et al Ill982h showed that natural 
thermal convection affects the mechanisms of momentum and heat transport from the wall and 
tends to flatten the velocity profile in t he bulk; these observatio n were confirmed by large-eddy 
simulations of a variable density fluid ( Zainali & Lessani 2010h . although with non-Boussinesq 
effects, as, e.g. profile asymmetries emerging at large stratifications. 

In this work we perform direct numerical simulations based on the lattice Boltzmann (here¬ 
after LB) method of an unstably (thermally) stratified turbulent channel flow. Our focus is on 
identifying the effects of buoyancy on the channel flow structure by comparing profiles of mean 
fields and fluctuations over a wide range of parameters with a pure (unstratified) channel flow. 
Our results show a decreased fluid throughput due to a strongly flattened velocity profile, which 
could, however, be fitted by a log-law with coefficients depending on the input controlling param¬ 
eters (friction Reynolds number Rct and Rayleigh number Ra, defined below), as also derived 
theoretically. 


2. Numerical method and simulation details 

We simulated a fluid enclosed between two walls, kept at fixed temperatures Th at ?/ = 0 and 
Tc = Th — a at y = 2H, and driven by a constant pressure gradient p' along the streamwise 
direction x (see figure [T]). The equations of motion are the incompressible Navier-Stokes equa¬ 
tion for the fluid velocity field it(x, t) in the Boussinesq approximation (namely, the density is 
assumed to be constant p = pQ, but for the linearised buoyancy term) 

dtu + u ■ Vu = —VP -I- vAu — /3g{T — Tm) + F, (2.1) 

{Tm = [Th + Tc )/2 being the mean temperature) coupled with the advection-diffusion equation 
for the temperature held r(x, t) 

dtT + u-VT = aAT. (2.2) 


In the above equations, P is the fluctuating pressure field (rescaled by po), g = —gy the ac¬ 
celeration of gravity, F = p'x a constant acceleration due to the imposed pressure head p'; v, 
a and /3 are the kinematic viscosity, the thermal diffusivity and th ermal expansion coefficient. 


respectively. As a numerical scheme, we adopted a 3d L B algorithm (iBenzi, Succi & Vergassola 


1992t Chen & Doolen 1998 ; Aidun & Clausen I 2()l()h with two probability d ensities (for den 
sity/momentum and for temperature, respectively) dHe. Chen & Doolenlll998 h. The method has 


been extensively used to study both th ermal conyection (IBenzi. Toschi & Trip iccione||1998 


Calzayarini, Toschi & Tripiccione 12002 ) and turbulent channel flow (Toschi et al\l 999t Toschi. Leyeaue & Ruiz-Chayarrfa 


2000t Biferale et al 2002 ): in particular iLayezzo, Clercx & Toschil (1201 Ih yalidated the code 
and tested grid resolutions against preyious studies with different numerical methods. 

We used a computational grid L x 2H x VL of 256 x 128 x 128 lattice points, with each run 
longer than 3 x 10® time steps (in LB units), in such a way to achieye statistically steady states 
of ^ 20QTl (Tl being the large-scale eddy turnoyer time). The domain under consideration is 
sketched in flgure[T] where a snapshot of an instantaneous held is shown. 

From equations (12.11) and (I2.21 i two dimensionless groups can be identified giying rise to two 
parameters which control the dynamics, namely the shear (or friction) Reynolds number, 

UtH 


RCr = 

V 

(Ur = Vp'H being the friction yelocity) and the Rayleigh number 
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Figure 1. Snapshots of the streamwise velocity field in the simulation domain (color map: velocity in the 
range [0; 20ut] from blue to red) for Rct = 205 and Ra = 0 (left) and Ra = 6.5 x 10® (right). The 
physical domain is also displaced along with the coordinate system adopted in the text. 
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Figure 2. Mean profiles of the streamwise velocity U{y), normalized with the friction velocity u-r, for 
various Rayleigh numbers (Ra = 0, Ra = 8.12 x 10®, Ra = 6.5 x 10® and Ra = 1.3 x 10^) at fixed 
shear Reynolds number Rcr- Left: Re^ = 130, Right: Rcr = 205. 


quantifying, respectively, the strength of the pressure induced shear and of the buoyancy with 
respect to viscous dissipation. 

We performed turbulent channel flow simulations with Rcr € [46,205]; for each Rcr we 
tuned the gravity (at fixed temperature jump), whence the buoyancy, spanning the range Ra G 
[0,1.3 x 10^]. 


3. Results 

Mean profiles of the streamwise velocity (7(y) = Thi (the overline indicates here and hence¬ 
forth averaging over planes parallel to the walls (•) = {LxLz)~^ f {■)dxdz) in the channel are 
shown in Figure ID For each Reynolds number {Rct = 130 and Rct = 205) we show data 
from simulations at three different Rayleigh numbers (Ra = 8.12 x 10®, Ra = 6.5 x 10® 
and Ra = 1.3 x 10^), besides the unstratifled channel flow (Ra = 0). An evident effect 
of thermal stratification is a decrease of the centreline velocity and a flattening of the pro¬ 
files at increasing Ra', indeed, mixing between the bulk and the boundary layer regions due 
to wall-normal thermal fluctuation s (or “plumes”) results in an increase of the effective wall drag 
( Hattori. Morita & Nagano Il2006h . Such effects are more pronounced for lower Rct, as it clearly 
appears by comparison of left and right panels of figure |2] These observations will be discussed 
more quantitatively under the light of the modelling in the next section. 

Figure[3shows the mean temperature profile with respect to y, for the same cases displayed in 
Figure |2] Here, we note the bending of the thermal profile in the bulk with increased Reynolds 
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Re^=130 Re,=205 



Figure 3. Mean profiles of the normalized temperature Q — 2{T — Tm)/A as a function of the height y 
for various Rayleigh numbers (Ra = 0, Ra = 8.12 x 10®, Ra = 6.5 x 10® and Ra = 1.3 x 10^) at fixed 
shear Reynolds number Rcr- Left: Rbt = 130, Right: Rct = 205. 
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Figure 4. Profiles of longitudinal — u% and transverse = Uy mean squared components of the 
fluctuating velocity field, normalized by the friction velocity Ut for various Rayleigh numbers (Ra = 0, 
Ra = 8.12 X 10®, Ra = 6.5 x 10® and Ra = 1.3 x 10^) at fixed shear Reynolds number Re-r- 


number, in contrast to the thermal shortcut observed in pure Rayleigh-Benard (at sufficiently high 
Rayleigh numbers) ( Scagliarini. Gvlfason & Toschi 12014 ). owing to the destruction or sweeping 
of the coherent plumes rising from the wall surfaces. 

Longitudinal and transverse mean squared components of the fluctuating velocity field (here¬ 
after Ui = Ui — Ui) are shown in Figure|4]for the same cases as in Figure |2] (normalized with the 
friction velocity). When the fluctuating quantities are concerned, most notably the lateral com¬ 
ponent is increased in the bulk region, and as the Rayleigh number is increased, the magnitude 
becomes comparable or greater than the stream-wise component. The near wall regions are also 
affected, primarily in the streamwise component, where an increase and a subsequent decrease 
in magnitude is observed as the Rayleigh number is increased. 

In order to quantify the overall effect of the thermal forcing on the channel flow, figure |5] 
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Figure 5. The skin-friction coefficient c/ = rm/(l/2[/o) as a function of the shear Reynolds number 7?er 
for various Rayleigh numbers Ra. The dashed lines are predictions obtained by approximating the mean 
velocity by a log-law over the whole channel; for the stratified cases we made use of equations 
with the same fitting parameter Cs = 2.5 for all Ra (see section|4]) 


shows the skin-friction coefficient c/ = Tu,/(1/2[/q) as a function of shear Reynolds number 
Rcr- The increased drag due to the thermal field is evident, resulting in higher than usual fric¬ 
tion coefficient. The effect is reduced asymptotically, as the Reynolds number is increased for a 
given Rayleigh number, emphasizing that shear becomes the dominant source of turbulence at 
sufficiently high Reynolds numbers. 


4. Modification of the law of the wall by buoyancy 

Inspection of the profiles of the mean quantities and the fluctuations showed that buoyancy 
alters the channel flow structure. We now focus on the mean velocity profiles. In Figure |6] we 
report the lin-log plot of profiles, in wall units, for Rcr = 185 and for various Ra. Convective 
motion, realised in the form of thermal plumes rising from the walls to the bulk, disturbs the 
coherence of the channel flow, resulting in a source of drag which decreases the mean velocity 
in the channel and hence the mass throughput, as well as flattening the velocity profile. 

In order to quantify these observations, in the following we attempt to generalize von Karman’s 
law of the wall accounting for buoyancy effects. To this aim, we propose a simple model con¬ 
structed from conservation laws and phenomenological arguments, which predicts that the vis¬ 
cous buffer is insensitive to buoyancy (i.e. [/+ = y~^), while in the log-law region the following 
relation 

U+{y+) = . )+B, (4.1) 

\l + Kcy+J 

holds (quantities have been expressed in wall units [/+ = U/ur and y~^ = y/S, with S = v/ur). 
In equation (14. Il l, k is the von Karman constant, B is an integration constant and nc is given by 

Rd 

Kc{Ra,Rer) = 

Cs is a phenomenological parameter. All profiles in Figure |6] agree with our predictions: all 
curves collapse on the function [/+ = j/+ for j/'*' < 10 (dashed lines) while in the logarithmic 
layer they can be fitted with equation (14.11) (dash-dotted lines), whose derivation is detailed in 
what follows. 
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Figure 6. Lin-log plot of the velocity profiles in wall units, for shear Reynolds number Rct = 185 (left 
panel) and Rayleigh numbers Ra = 0 (^=), Ra = 8.12 x 10® (•) and Ra = 6.5 x 10® (■). The dashed 
lines indicate the linear law, valid in the viscous layer, — y^, while the dash-dotted ones correspond to 
the prediction | |4.U -( I4.2> with k = 0.42, Cs = 2.5 and B € [5.8, 6.7] (Pr = 1). 


Applying the Reynolds decomposition u(x, f) = U{y)x + u(x,f) for the velocity field and 
averaging the Navie r-Stokes equ ation for the x velocity component, we get, upon integration in 
y, an exact relation (IPope 1120011) between the mean shear S{y) = dU{y)/dy and the Reynolds 


stress Txy{y) = —UxUy namely: 


vS{y) + Txy{y) =p'{H -y), 


or, in wall units: 


5+(y+) 


■xy(y^) (1 


(4.3) 


(4.4) 


Notice that equations (I4.3l )-( l4.4b remain valid also in the presence of buoyancy, since the forcing 
term /3p(T — Tm)y does not act along the streamwise direction. For very large Rcr (in principle 
Rct oo), equation (l4~4l i gives + T^y = 1; close to the wall (for ?/+ < 10), then, where 
viscous terms dominate, S'^ <C r^y whence S'^ ~ 1 which implies [/+ ~ y~^, as anticipated 

above. Conversely, away from the wall (ty pically for y^ ^ 3 0), turbulent fluctuations dominate 
over viscous processes and T+ <C 5'+, i.e. dUvov et a! II2004I) 


•xy 


const: 


(4.5) 


such result alone, however, does not give any further insight on the behaviour of the velocity 
profile in the corresponding region of the channel. To th is aim, we n eed to consider the budget 
equation of turbulent kinetic energy Ex{y) = |up/2 dPope 2001 ). In the log-law layer the 
energy prod uction V = Txy{y)S{y) + I3g(j){y) is balanced by dissipation e{y) = 2v{diUjY 
dPope II2OOII) : the latter cannot be calculated exactly but only inferred by means of phenomeno¬ 
logical arguments: for large Reynolds number and outside the viscous boundary layer, the dissi¬ 
pation equals the turbulent energy flux, which can be estimated as E kIv)/ riy), where T{y) oc 
y / \/EK{y) is the typical eddy turn-over time at y ( L’vov ef^l2004l) . Thus, the energy balance 
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Figure 7. Main panel: Ratio of the turbulent kinetic energy Ek over the Reynolds stress Txy'. the plateau 
attained for t/^ ^ 40 suggests that the two quantities are proportional to each other in the log-law region. 
Inset: Ratio of the turbulent heat flux over the mean shear; the solid line indicates a linear dependence on 
and represents a numerical check of the validity of the assumptions leading to equation d4.1 lb . 


equation reads: 

p3/2 / \ 

° ^ ^ = Txy{y)S{y) + figfiiy), (4.6) 

y 

where a is a non-dimensional parameter of order unity. 

Equation (14.6b is not yet closed: a relation between Ek {y) and r^y (y) is required. Dimensional 
analysis suggests that the two quantities should be proportional to each other, i.e. EK{y) oc 
Txy(.y) (IL’vov et aZll2004l) : such assumption can be readily tested in the numerical simulations: 
in the main panel of Figure|7]we show that the ratio Ek {y) j^xy {y), indeed, approaches a constant 
value for ~ 40. Equation (14.6b can be then rewritten as (in view also of (14.5b ) 

— = cEKS{y) + figfiiy), (4.7) 


with c yet another dimensionless number. The energy production by buoyancy te rm, viz, the heat 


flux, needs also to be modelled; althoug h more refined closures can be employed ((.Tohansson & Wikstrom 
1999t Hattori. Morita & Nagano 20061) involving tensorial eddy diffusivities and coupling with 


gradients of the temperature field in all directions, we adopt a simple standard mixing length 
ansatz, i.e. 


4>{y) = UyO = -il^dyUdyT. (4.8) 

A first order closure like (14.8b must not be expected to work well for second order quantities like 
temperature fluctuations, Reynolds stresses, etc, but, as we will show, mean streamwise velocity 
profiles are satisfactorily reproduced through such model; according to Prandtl’s hypothesis the 
mixing length Era is proportional to the distance from the wall, i.e. Im oc y, whence 


(j){y) = UyO = -by'^dyUdyT, 


(4.9) 


where 6 is a numerical constant. The mean temperature gradient will, in principle, depend it¬ 
self on the shear profile; however, in a perturbative spirit, we postulate here, for simplicity, a 
logarithmic form, such that 


dT 


dy 


h 

y 


(4.10) 
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Figure 8. Plots of the parameter K.c{Ra, Rct) of the model as function of Ra for fixed Rct (main panel) 
and as function of Rbt for fixed Ra (inset). The dashed lines represent equation O with constant Cs 
fixed at Cs = 2.5 for all cases. 
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where T* = 


:)tA 


2H Ut I * I I I - 

lates between the p assive scalar case Llohansson & Wikstrom I (119991) and the natural convection 


is a scale for the dynamic temperature. The for m (14.10b somehow interpo- 


Ahlers et al (120121) : inserted in the expression for the heat flux (14.9b it gives 


rp{C) 

4>{y) cx by^Siy)^ = h'yS{y)TP. (4.11) 

y 

For the sake of validation of our arguments, we check equation (14.1 lb (which predicts (/)(?/) oc 
yS{y)) against the numerics in the inset of Figure |7] finding a reasonably good agreement. In¬ 
serting (14.1 lb into (14.7b provides 


y 


cEKS{y) + /^9b'yS{y)^j^ 


which can be recast, introducing the wall units and the definitions of Rct and Ra, in the following 
form 


1 

Ky+ 


S+ 


l + Cs 


Pr^Re^ ) ’ 


(4.12) 


it must be noticed that the two parameters k and Cs are combinations of dimensionless quantities 
emerging in the derivation, which cannot be, however, derived from first principles; fits of the 
numerical data indicate that good estimates are the values k = 0.42 and Cs = 2.5. Explicitating 
the shear term in (14.12b . we get 


dU+ _ 1 

dy+ Ky+{1 + Kcy~'')’ 


(4.13) 


whose integration Anally yields expressions ( 14.1b and ( 14.2b for the velocity profile. The robust¬ 
ness of the model is confirmed in Figure|^where we plot the fitted values of the parameter kq as 
function of Ra for fixed Rer (main panel) and as function of Re^ for fixed Ra (inset), together 
with the predictions of equation (14.2b (dashed lines). 


5. Conclusions 

We have studied, by means of direct numerical simulations based on a thermal lattice Boltz¬ 
mann algorithm, the dynamics of a turbulent channel flow under a gravity field orthogonal to 
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the streamwise direction coupled to an imposed temperature difference between the top (cold) 
wall and the bottom (hot) wall. The resulting unstably stratihed conhguration flattens the ve¬ 
locity prohle and decreases the centreline value when the buoyancy strength is increased. This 
effective drag shows up in an enhancement of the friction coefficient. The action of buoyancy on 
the boundary layer structure has also been probed looking at other relevant statistical quantities 
in wall bounded turbulent system, such as Reynolds stress; we have found that, as the Rayleigh 
number is increased, the squared wall normal velocity grows in the bulk becoming comparable 
or even larger than the streamwise component (which, in turn, is depleted close to the wall). To 
provide a quantitative interpretation of the numerical hndings, we have proposed a phenomeno¬ 
logical model resulting in a modihed logarithmic law of the boundary layer; such model could 
successfully capture the various velocity prohles at changing the shear Reynolds and Rayleigh 
number, with just one adjustable parameter. 
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